
%%%%%%%%%%%%%%%%%
% Factor model for options market making
% Liuren Wu, liuren.wu@baruch.cuny.edu
%%%%%%%%%%%%%%%%%
nx=6; nxe=nx+2; nlag=12;
minN=200;winth=0.01; 

dds=Z.dds;
ids=Z.ids;
dd=Z.dd; T=length(dd);
IV=mean(Z.IVV(:,1:2,1),2); %ATMV
y=mean(Z.NRV(:,1:2,1,3),2)*100; %straddle return, daily updated hedged
OISV=Z.ControlVariables(:,5);% Option open interested scaled by stock volume
x=[Z.Factors, OISV]; %Factors from return factor model (HC, VR, JR) + Option open interest
bas=Z.OBASA*100; %bid-ask spread (average of call and put)
HVV=Z.HVV; %historical vol at 3 horizons



Bfv=NaN(T,nx);Afv=NaN(T,1);
statv=NaN(10,nx+2,T);
R2v=NaN(T,1);
VB=NaN(T,nx);
parfor t=nlag+1:T
    [Af,Bf,R2s,stats,vbt] =FunDailyORSestimationBAS(t,y,x,bas,HVV,IV,ids,dds,dd,nlag,nx,winth,minN);
    
    Bfv(t,:)=Bf; R2v(t)=R2s;
    Afv(t)=Af;
    statv(:,:,t)=stats;
    VB(t,:)=vbt;
    
end

fl=[Afv,Bfv,100*R2v];
Tf=sum(isfinite(fl));
flstats=[nanmean(fl)', nanstd(fl,0)',sqrt(Tf)'.*nanmean(fl)'./nanneweystd(fl)',nanacf(fl,1)', skewness(fl,0)', kurtosis(fl,0)'-3] ;



xnamel={'Intercept', 'Hedging cost','Volatility risk','Jump risk','Open Interest', 'Historical premium','Volatility premium','Adjusted $R^2$'}';

ax=reshape(statv(1,1,:),T,1); Tf=sum(isfinite(ax));
axs=[nanmean(ax), nanstd(ax,0)',sqrt(Tf)'.*nanmean(ax)'./nanneweystd(ax)', nanacf(ax,1)', skewness(ax,0)', kurtosis(ax,0)'-3] ;

fprintf(1, '  \n');
fprintf(1, '%40s ', 'Average');fprintf(1, [repmat(' & %8.2f  ',1,2),repmat(' & %8.2f  ',1,4),' \\\\ \n'], axs);

mVB=[NaN;100*nanmean(VB)'; NaN];


fprintf(1, '\\\\  \n');
for j=1:1;fprintf(1, '%40s ', xnamel{j});fprintf(1, [repmat(' & %8.2f  ',1,2),repmat(' & %8.2f  ',1,4), ' \\\\ \n'], [flstats(j,:)]);end
for j=[2,3,4,6,7];fprintf(1, '%40s ', xnamel{j});fprintf(1, [repmat(' & %8.2f  ',1,2),repmat(' & %8.2f  ',1,5), ' \\\\ \n'], [flstats(j,:),mVB(j)]);end
fprintf(1, '\\\\  \n');
for j=[5];fprintf(1, '%40s ', xnamel{j});fprintf(1, [repmat(' & %8.2f  ',1,2),repmat(' & %8.2f  ',1,5), ' \\\\ \n'], [flstats(j,:),mVB(j)]);end
fprintf(1, '\\\\  \n');
for j=8;fprintf(1, '%40s ', xnamel{j});fprintf(1, [repmat(' & %8.2f  ',1,2),repmat(' & %8.2f  ',1,4), ' \\\\ \n'], [flstats(j,:)]);end

